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NATIONAL  ADVISORY  COMMITTEE  FOR  AERONAUTICS 


TECHNICAL  MEMORANDUM  1^03 


ON  THE  INSTABILITY  OF  METHODS  FOR  THE  INTEGRATION 

OF  ORDINARY  DIFFERENTIAL  EQUATIONS1 

By  Heinz  Rutishauser^ 


In  spite  of  the  remarkable  publication  of  J.  Todd  (ref .  l)  the 
essential  points  of  which  are  related  below,  the  author  has  since 
observed  several  times  methods  for  the  nummerical  integration  of  differ- 
ential equations  which,  although  subject  only  to  a  temptingly  small 
truncation  error,  nevertheless  involve  the  great  danger  of  numerical 
instability.   I  want  to  state  beforehand  that  this  danger  hardly  exists 
for  the  well-tested  methods  of  Runge-Kutta  and  Adams  (extrapolation 
methods)  if  they  are  applied  correctly. 

It  is  a  natural  characteristic  that  a  differential  equation  to  be 
solved  numerically  is  approximated  by  a  difference  equation,  and  that 
the  latter  is  then  solved.   In  order  not  to  be  forced  to  select  an  all 
too  small  interval,  one  prefers  difference  equations  which  approximate 
the  differential  equation  as  closely  as  possible  but  in  compensation 
are  of  higher  order  than  the  original  differential  equation.  Precisely 
in  this,  however,  there  lies  a  danger  because  the  difference  equation 
thereby  has  a  greater  diversity  of  possible  solutions,  and  it  may  well 
happen  that  the  numerical  integration  yields  precisely  one  of  the 
extraneous  solutions  which  only  at  the  beginning  is  in  any  way  related 
to  the  desired  solution  of  the  differential  equation.   In  the  paper  of 
J.  Todd  mentioned  before  several  examples  of  this  type  have  been 
enumerated . 

Consideration  of  the  pertinent  variation  equation  is  particularly 
informative.   It  is  very  well  possible  that  the  differential  variation 
equation  is  stable,  that  is,  that  it  contains  only  converging  solutions, 
whereas  the  difference  variation  equation  is  unstable  since  it  possesses, 
due  to  the  increased  diversity  of  solutions,  aside  from  the  converging 
solutions,  also  solutions  which  increase  exponentially.  A  deviation 
from  the  correct  solution,  once  it  exists,  small  as  it  may  be  -  and  such 

Uber  die  Instabilitat  von  Methoden  zur  Integration  gewohnlicher 
Dif f erentialgleichungen , "  ZAMP,  Kurze  Mitteilungen,  vol.  Ill,  1952, 
pp.  65-7^. 

^Institut  fur  angewandte  Mathematik  der  ETH,  Zurich. 
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deviations  are  unavoidable,  because  of  the  rounding-off  errors  -  there- 
fore increases  rapidly  and  may  finally  falsify  considerably  the  solution 
obtained.  Yet  -  we  want  to  emphasize  this  once  more  -  this  instability 
is  caused  only  by  an  inappropriate  integration  method. 

In  the  following  discussion,  several  customary  methods  are  examined 
from  this  viewpoint,  and  simple  criteria  for  the  stability  of  such  meth- 
ods are  indicated.  For  the  rest,  this  report  does  not  deal  with  error 
estimates. 


DIFFERENTIAL  EQUATIONS  OF  THE  FIRST  ORDER 
y'  =  f(x,y) 


Variation  equation 


T)1  =  M  TJ 


(a)   Integration  by  Means  of  Simpson's  Rule-3 

yk+i  -  yk-i + 1  (y'k+i +  k*\ +  y'k-i)  (1) 

This  relationship,  together  with  the  differential  equation,  yields 
two  equations  for  the  unknown  quantities  Jv,-,      and  y'     which  are 

solved  mostly  by  iteration.   If  the  differential  equation  is  linear  or 
quadratic  in  y,  the  iteration  can  be  avoided.  We  assume,  however,  that 
one  passes  over,  in  any  case,  to  the  next  integration  interval  only  when 
the  relationship  (l)  is  satisfied. 

The  difference  variation  equation  pertaining  to  (l)  is  evidently 


The  phenomenon  was  observed  on  this  example,  and  correctly  inter- 
preted, also  by  Mr.  G.  Dahlquist,  Stockholm  (Lecture  at  the  GaMM- 
convention  1951  &t  Freiburg  im  Breisgau) .   Compare  also:   Z.  angew. 
Math.  Mech.,  31,  239,  1951- 

^Where  h  signifies  the  length  of  the  integration  interval,  and 
y.   stands  as  an  abbreviation  for  y(kh). 
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1*1(1  - 1  fy,k+1)  -  f  fy,k\  -(i  ♦  |  vjvi  -  ° 

If  one  assumes  f   to  be  constant,  and  chooses  the  expression  t],  =  A 

for  the  solution  of  this  equation,  one  obtains  for  A  a  quadratic 
equation  with  the  solutions 

h2   2  hf„ 


k 


\  =   1  +  hf  +  —  fy  +  .  .  .  «,   e  y 

One  recognizes  easily  that,  of  the  two  fundamental  solutions  T),  ,  =  A]_ 
and  ^2  k  =  ^2   °^  the  difference  variation  equation,  the  first  one 

approximates  the  solution  of  the  differential  variation  equation  whereas 
the  second  one  is  brought  in  by  the  numerical  method. 

In  particular 

A2  ~  (-1)  e    Ji 

represents  for  fy<0,  thus  precisely  when  the  differential  equation  is 

stable,  an  oscillation  which  is  slowly  exponentially  increasing.  This 
has  the  effect  that  a  small  disturbance  of  the  numerical  solution  - 
caused  by  a  rounding-off  or  a  truncation  error  -  is  intensified  in  the 
further  course  of  the  integration  and  finally  gets  completely  out  of 
hand.   In  Collatz1  (ref.  2)  book,  the  phenomenon  is  denoted  as  "roughening 
phenomenon";  means  for  elimination  of  this  inconvenience  are  given.  On 
the  other  hand,  the  explanation  given  there  is  not  complete;  the  phe- 
nomenon occurs  only  for  fy  <  0;  there  is  nothing  to  be  apprehensive  of 

for  f y  =  0  which  is  very  important  particularly  in  regard  to  the  ordi- 
nary Simpson's  integration  rule  (fy  =  0). 

(b)  Integration  According  to  Runge-Kutta  and  Similar  Methods 

Since  these  methods  calculate  y,  ■,   from  y^  according  to  a 
prescribed  rule  and  without  use  of  the  preceding  values  Y\_\>   ^-2'  •  *  •> 
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the  order  remains  unchanged  in  the  transition  from  the  differential 
equation  to  the  difference  equation;  thus  no  foreign  solutions  are  brought 
in,  and  no  instability  is  to  be  feared. 

The  same  property  can  be  found  in  a  method  indicated  by  W.  E.  Milne 
(ref.  3). 

(c)   Integration  According  to  Adams 
We  consider  a  four-point  formula 

yk+i  =  yk  +  ^  fer'n+i  +  19^'k  -  5y'k_i  +  y\_2)  +  n?  .  .  .       (2) 

This  again  yields,  together  with  the  differential  equation,  two  equa- 
tions for  the  unknown  quantities  y.  ■,   and  y'k+-i   which  are  mostly 
solved  by  iteration. 

The  difference  variation  equation  pertaining  to  (2)  becomes 


^  -  f  fy,fcn)%n  -  (x  +  lr  fy,k)nk  +  S  ^k-i^-i  -  ^  c  v_^>_5  =  0 


^  1y,k-2T)k-2 


If  one  again  considers  f„  as  constant,  the  expression  T],  =  A   yields 

an  equation  of  the  third  order  for  A.  A  solution  A,   of  this  equation 

hfy  v 

lies  very  close  to  e   ,  therefore  ^l  k  =  ^1   corresponds  to  the 

'  k  k 

solution  of  the  differential  variation  equation  whereas  Ag   and  A* 

are  extraneous  solutions . 

However,  the  equation  for  A  is  reduced  to  A"*  -  A^  =  0  when  h 
tends  toward  0  so  that,  for  a  sufficiently  small  h,  one  will  have  at 


any  rate  small  A2  and  A*,  namely  ~  ±\|-hfy/21+.  The  extraneous  solu- 

k  k 

tions  T]2  k  =  ^2   and-  Ix  k  =  'N   thus  converge  rapidly.  For  a  suffi- 
ciently small  h,  Adams'  method  is  therefore  stable. 

(d)  Variants  of  (c) 

In  order  to  improve  the  accuracy  of  Adams '  method,  one  may  use  also 
other  expressions  for  the  corrector  instead  of  (2).  As  long  as  the 


•   •   • 
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corresponding  five-  or  six-point  formulas  are  involved,  there  are  no 
objections,  but  one  has  to  be  careful  when  yk+l  is  not  calculated 

from  yk  and  the  derivatives  as  in  (2)  but  perhaps  from  yv_]_  or 
yk-3  an(^  ^e  derivatives,  as  for  instance  in 

yk+i  =  yk-3  +  ff  (7y'i-+i  +  52y'k  +  I2y'k_!  +  32y'k_2  +  ly\.3\  +  hT 

In  fact  the  pertinent  difference  variation  equation  has  a  solution 

Ifc.  Ak   With   X=-(l-^fy+.   •  >j 

so  that  the  method  is  unstable  for  fv  <  0. 

DIFFERENTIAL  EQUATIONS  OF  THE  SECOND  ORDER 

y"  =  f(x,y,y') 

Insofar  as  these  equations  are  solved  by  separation  into  a  system 
of  two  equations  of  the  first  order,  what  was  said  so  far  is  valid. 
Particularly  in  the  case  of  numerical  integration  of  damped  oscillations 
we  must  caution  against  the  methods  (a)  and  (d) . 

However,  there  exist  also  methods  which  solve  an  equation  of  the 
second  order  without  transformation  into  a  system: 

(e)  The  Method  of  Central  Differences5 

The  formulas  on  which  this  method  is  based  are  (especially  for 
second  order) 


yk+i  =  2yk  -  y-k-i  +  ^(y'W  +  10A  +  A-i) 
yk+i  =  y,k-i +  !(y,,k+i +  ^\  +  y"k-i) 


(5) 


^Compare  reference  2,   p.  80. 
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They  yield,  together  with  the  differential  equation,  three  equations 
for  the  unknown  quantities  yv+]_>   y'k+i>  an<i  y"k+l"  ^he  two  simul~ 

taneous  difference  variation  equations  pertinent  to  (k)    and  (5)  are 
solved  with  the  expression  tj,  =  pA  ,  r\'      =  qA  ,  under  the  assumption 

of  a  constant  fv  and  f ,}   because  of 
j  y 


'Vl-  fyVi+  V^'k+i 

one  obtains,  with  the  abbreviations  a  for  h2f v/l2  and  b  for 
hf  1/3,  the  equations  ' 

1/ 


pKA2   -  2A  +  l)    -  a(A2  +  10A  +  l)    - 
q  -  b(A2  +  10A  +  1)   =  0  [from  (k)\ 


(A2  -  1)  -  b(A^  +  i+A  +  1) 


p  £  a(A2  +  k~h   +  1)  =  0     [from  (5) 


(6) 


These  equations  can  exist  simultaneously  with  (p,q)  ^  (0,0)   only  when 
the  determinant  of  this  equation  system  for  p  and  q  vanishes;  one 
obtains  after  some  calculations 


(A  -  1)  ,(A2  -  1)(A  -  1)  -  a(A  +  1)(A2  +  10A  +  l) 


)(A  -  1)(A2  +  1+A  +  1)]  = 


0 


The  four  solutions  of  these  equations  are 


A]_  =  1  +  a]  h  +  .  . 

A2  =  1  +  a^h  +  .  . 

> 
A3  =  1 

\   =  -fl-|fy.  +  .  .  ) 


where  a-j_  and  a,o  are  the  solutions  of 
the  equation  a2  -  afv>  -  f 


y 
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Evidently  A]_k  and  Apk  are  the  regular  solutions  of  the  difference 

variation  equation;  they  correspond  to  two  fundamental  solutions  of 
the  differential  variation  equation;   A*  and  A^,  in  contrast,  are 

extraneous.  As  long  as  f  1  ^  0,  there  is  nothing  to  fear,  in  partic- 
ular, the  method  may  be  strongly  recommended  for  a  y'-free  equation, 

ed  oscillations) 

%k  =  A4k.  (-l)V(W3)fy. 


but  for  fyi  <  0  (damped  oscillations) 


increases,  and  A*,  too,  may  still  become  dangerous  because  A*  ^  =   1 

also  becomes  finally  very  large,  compared  to  a  function  converging  to- 
ward zero. 

The  author  completely  calculated  the  example  y"  +  y1  +  1.25y  =  0 
with  the  initial  conditions  y  =  0,   y'  =  1  (exact  solution: 

¥      \  ° 

e  sin  x/  on  the  sequence-controlled  computing  machine  of  the  ETH. 

There  follow  a  few  excerpts  from  the  thus  obtained  table  of  func- 
tions (we  calculated  with  h  =  O.l): 


X 

y 

y' 

k.8 

k.9 

5-0 

5.1 
5.2 

-0.0905699 

-  .0847792 

-  .0787152 

-  .0722891 

-  .0656175 

0.0551227 
.0581+842 
.06261+10 
.0656575 
.0676070 

In  this  region  nothing  conspicuous  is  noticeable  yet,  the  y-values 
deviate  from  the  true  values  approximately  by  one  in  the  last  decimal 
place,  and  only  formation  of  the  differences  for  the  y' -values  reveals 
a  certain  irregularity.  For  t  =  17,  however,  the  influence  of  A^ 

becomes  pronounced  for  the  y' -values  and  also  for  the  differences  of 
the  y-values : 


8 
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X 

y 

y" 

17.0 

-0.00019574 

0.00005017 

17.1 

-  .00019061 

.00005253 

17.2 

-  .00018366 

.00008620 

17.3 

-  .0001752^ 

.00008239 

17.4 

-  .00016548 

.00011235 

17.5 

-  .00015475 

.00010258 

The  considerably  weaker  oscillation  of  the  y-values  follows  also 
from  the  equations  (6):   for  A  =  \      one  obtains  from  the  first  of 
these  equations 


—  ~  —  D 


8      h 
'—  =  -  —  f,_|.  here  therefore  p  ~  — *- 

4  6     y  '  600 


The  further  course  of  numerical  integration  does  not  require  any  comment: 


X 

y 

y' 

22.8 

-0.00000815 

0.00005320 

22.9 

-  .00000864 

- .00006078 

23.0 

-  .00000862 

.00005968 

23.1 

-  .00000887 

-.00006247 

23-2 

-  .00000861 

.00006601 

23-3 

-  .00000868 

- .00006486 

29.5 

-  .00000140 

- .00053037 

29.6 

.00000041 

.00054889 

29.7 

-  .00000144 

- .00056693 

29.8 

.00000050 

.00058682 

29-9 

-  .00000148 

.00060603 

30.0 

.00000060 

- .00062735 

The  author  is  well  aware  that  the  assumption  of  a  constant  fy  and 
fvi   in  the  above  considerations  greatly  restricts  the  generality.  How- 
ever, the  results  show  that  what  matters  is  only  the  sign  of  these  quan- 
tities, and  this  sign  is  indeed  invariable  in  a  great  many  cases.  The 
statements  are,  therefore,  qualitatively  almost  generally  valid.   Only 
when  fy  changes  its  sign,  from  time  to  time  in  the  course  of  the  inte- 
gration, for  instance  when  method  (a)  is  being  used,  a  special  case 
arises  since  the  occurring  oscillations  alternately  increase  and  are 
damped  again. 
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GENERAL  CONSIDERATIONS 


Consideration  of  the  examples  suggests  the  conjecture  that  insta- 
bility may  occur  precisely  in  the  case  of  integration  methods  which 
form  y^+x  by  integration  of  y1   over  several  intervals  (two  intervals 

in  method  (a),  one  interval  in  (b)  and  (c),  four  intervals  in  (d),  two 
intervals  in  (e)).  However,  this  is  not  exactly  the  case  and  we  shall 
therefore  subject  a  general  integration  method  to  an  examination. 

Almost  all  known  methods  use  relationships  which  are  contained  in 
the  general  formula 


'k+l 


0 

H 


Oj   •'k+j  Z_     lj  J    k+j 

-m  -m 


(0)    (N) 


h      >     a,  .   y,  ,  . 

-m 


k+l  "  Z_  alj  y  k+j  + 

-m  -m 


h  £.  4JV\+J  ♦ . . .  + 


**  t « 


-m 


-m 


N-n+1 


-m 


nj 


k+j 


(n-1)    (N) 


^j 


-m 


k+j 


; 


(7) 


n-l=v^a(n-l)    (n-1)   +  h  ^  a(n-l)    (n) 
k+l       Z_     n-l,,fk+i  /=_     n.i       Jk+j 


10 
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There  are,  in  addition,  differential  equations 


■n  k+i^k+i'  •  •  •>  yk+i  =  ° 


n+1 


(xk+i>  • 


k+i'    * 


in+1)' 


?w 


\xk+l'yk+l' 


'    •'   yk+l 


=  0 


(N) 
rk+l 


)• 


0 


(8) 


which  form,  together  with  the  equations  (7),   N+1  equations  for  the 

(N) 
N  +  1  unknowns  Yv+i>  ^  k+l>  *  "  •>   V     ■<•      Generally,   N  =  n;  however, 

W.  E.  Milne  (ref .  3)  uses  in  the  method  previously  mentioned  higher 
derivatives  than  those  appearing  in  the  differential  equation.  There- 
fore he  differentiates  them  several  times  in  order  to  obtain  the  required 
number  of  relationships.  One  may  thus  obtain  the  equations  F  ,  .  .  . 

to  F^j  by  differentiating  the  initial  equation  Fn. 


If  one,  furthermore,  brings  everything  in  the  equations  (7)  to  one 
side,  the  variation  equations  for  the  entire  system  (7)  and  (8)  read 


evidently  I with  a.  ,  =  -1] 


a(i)hH-i>)  =  0 


H=i  j=-m 


HJ 


k+j 


(i  =  0,  1,  .  .  .,  n  -  l) 


^±      >)   . 

T],     =0 


0   3y(n)    X 


(i  =  n,  n  +  1,  .  .  .,  N) 


(n)     i 
If  one  uses  for  this  the  expression  T) .   =  p  V,  one  obtains  in  sub- 
stituting a  system  of  N+1  simultaneous  equations  for  the  p   which 


TM  1U03 


11 


can  be  satisfied  only  when  the  determinant  of  the  system  disappears 


|i=i\j=-m  ^°   / 


=  0      (l  =  0,  1,  .  .  .,  n  -  1) 


X   bf 


0  oy(^ 


(i  =  n,  n  +  1,  .  .  .,  N) 


If  one  defines,  in  addition,  with  the  coefficients  a  .   appearing  in 
the  formulas  (7)  the  functions 


Aiu(A)  =  5   a(i)AJ+m 


wherein  A.  =  0  for  |i  <  i,  the  characteristic  determinant  reads  as 
follows 


D(A)  = 


AOO 
0 


hAd 
All 


h2Ao2 


hA 


12 


0        0 


0 


An-l^n-l 


hNAQN 


hN-n+lA  ,  „ 
n     An-1,N 


*n,y     *n,y' 
Fn+l,y   Fn+l,y' 


Fn,y(n)   0       0 


FN,y' 


N,y 


\yW 


12 
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If  the  method  is  to  reproduce  exactly  a  polynomial  of  the  ith  degree, 
which  is  the  solution  of  a  differential  equation,  together  with  its 
derivatives  -  and  this  one  may  require  -  the  conditions 

y(l)=  1     y(i+l)  =  y(i+2)  =  y(i+3)  =  .  .  .  =  y(N)  =  0 

J  J        .3        J  j 

must  be  compatible.  Hence  follows,  however,  by  substitution  into  the 
ith  of  the  equations  (7):V~  a^1)  =  0,  therefore  A.,.(l)  =  0. 

The  equation  D(A)  =0  which  is  decisive  for  the  stability  of  the 
method  must  have  n  solutions  in  the  neighborhood  of  A  =  1,  corre- 
sponding to  the  n  independent  solutions  of  the  variation  equation. 
In  fact  one  finds  for  h  =  0  where  D(A)  is,  except  for  one  factor, 
reduced  to  AqcAu  "  '  '  ^n-1  n-1   "that  A  =  1  is  an  n-fold  zero  of 

D(A),  because  of  Ai;L(l)  =  0. 

All  other  zeros  of  D(A)   correspond  therefore  to  extraneous  solu- 
tions of  the  difference  equation;  in  order  to  make  the  method  stable, 
they  must  lie,  for  a  sufficiently  small  h,  in  the  interior  or  at  most 
on  the  periphery  of  the  unit  circle.  This  is  certainly  the  case  when 
for  h  =  0  all  zeros  of  D(A)  lie  in  the  interior  of  the  unit  circle, 
and  certainly  not  when  individual  ones  are  outside  it.  Therefore: 

Sufficient  condition  for  the  stability  of  the  method  (7)  for 
sufficiently  small  h: 

All  functions 

Aii(A)  =7"a(iKJ+m 

-m   ^ 

possess,  aside  from  the  trivial  simple  zero  A  =  1  only  zeros  with 

|a|   <  1. 

Necessary  condition:   None  of  the  functions  A^(A)   has  a  zero 
outside  the  unit  circle. 

<J.  Todd  considered  methods  which  satisfy  not  even  the  necessary 
condition.  The  solution  obtained  then  becomes  completely  useless 
already  after  a  few  intervals . 
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The  remaining  functions  A.  (A)  can  influence  the  stability  only 

when  the  necessary  condition,  but  not  the  sufficient  condition,  is 
satisfied. 


APPLICATIONS 


For  the  formula  (l)  there  results  (one  has  N  =  n  =  1,  m  =  l) 


Aqo  =  -A2  +  1      Aq!  =  J  (A2  +  k\  +   1)      F,  =  y«  -  f(x,y) 

3 


Therefore 


D(A)  = 


1  _  A2    £  (A2  +  4A  +  1) 
3 


-fu 


The  fact  that  Aqq  has  two  zeros  with   |a|  =  1  already  suggests  cau- 
tion, but  moreover  one  reads  off  imediately  that  D(A)  is  positive 
for  f y  <  0  and  A  =  -1,  and  negative,  in  contrast,  for  A  =  -<». 

Thus  a  zero  lies  to  the  left  of  -1;  the  method  is  unstable. 

For  the  formula  ^.k2   in  the  book  of  Collatz  mentioned  (p.  8l), 
there  is   (N  =  n  =  k,      m  =  l) 


?k  =   yIV  -  f(x,y,y',y",ym) 
Aqq  =  A22  =  -(A  -  l)2 

Ali  -  A33  =  -^  +  1 
A^  =  2A 


^ 


=  1  (A2  +  UA  +  1) 


All  other  A.   occur  only  with  at  least  h2  in  the  determinant, 
one  has,  except  for  terms  with  h^ 


Thus 


ih 


TM  1U03 


D(X)  = 


•(A  -  I)2 
0 
0 
0 


-A2  +  1 
0 
0 


-v 


0 

2Ah 

■(A  -  I)2 

0 


0 

0 

0 

-A2  +  1 


V 


0 

0 

0 

h  (A2  +  4A  +  l) 
3 


For  A  =  -oo,   D  is  positive,  for  A  =  -1  -  e,   D  has  the  sign  of 


-A2  +  1    I  (A2  +  k\  +   1) 
3 


y. 


Which  for  f m  <  0  and  a  sufficiently  small  e  is  negative .  Therefore 
this  method  is  unstable  for  f _jit  <  0. 


On  the  other  hand  it  is  easy  to  indicate  methods  which  are  always 
stable.  One  need  only  shape  the  formulas  (7)  in  such  a  manner  that 
every  line  begins  with 


.yfi)   =  yU)   +  h  Y~   a(i)   .y(i+l)   + 


yk+l  "  yk 


-m 


"i+l,jy  j+k 


(i  =  0,1,  .  .  .,  n-1) 


.m+1  ,  ,m 


Thereby  A. .(A)  =  -A    +  Am  and  has  therefore  only  the  trivial  zero 
A  =  1  on  the  periphery  of  the  unit  circle . 


SUMMARY 


In  the  numerical  solution  of  a  differential  equation  as  a  difference 
equation,  the  latter  is  usually  of  higher  order  and  therefore  has  more 
solutions  than  the  original  differential  equation.   It  may  well  be  tnat 
some  of  these  "extra"  solutions  grow  faster  than  any  solution  of  the 
given  equation;  in  this  case  the  computational  solution  has  the  tend- 
ency to  follow  one  of  these  and  has  after  a  certain  number  of  integra- 
tion steps  nothing  to  do  with  the  original  differential  equation. 
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The  author  gives  some  examples  and  a  criterion  for  stability  of 
integration  methods.  This  criterion  is  then  applied  to  some  well-known 
integration  formulas . 


Translated  by  Mary  L.  Mahler 
National  Advisory  Committee 
for  Aeronautics 
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